Efficient ab initio calculations of bound and continuum excitons 
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We present calculations of the absorption spectrum of semiconductors and insulators compar- 
ing various approaches: (i) the two-particle Bethe-Salpeter equation of Many-Body Perturbation 
Theory; (ii) time-dependent density-functional theory using a recently developed kernel that was 
derived from the Bethe-Salpeter equation; (iii) a scheme that we propose in the present work and 
that allows one to derive different parameter-free approximations to (ii) . We show that all methods 
reproduce the series of bound excitons in the gap of solid argon, as well as continuum excitons in 
semiconductors. This is even true for the simplest static approximation, which allows us to refor- 
mulate the equations in a way such that the scaling of the calculations with number of atoms equals 
the one of the Random Phase Approximation. 

PACS numbers: 71.10.-w, 78.20.Bh, 71.35.-y, 71.15.Qe 
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Time-dependent density-functional theory (TDDFT) 
[l[ is more and more considered to be a promising ap- 
proach for the calculation of neutral electronic excita- 
tions, even in extended systems 0, In linear response, 
spectra are described by the Kohn-Sham independent- 
particle polarizability Xo^"^ and the frequency-dependent 
exchange-correlation (xc) kernel fxc- The widely used 
adiabatic local-density approximation 0, [3 (TDLDA), 
with its static and short-ranged kernel, often yields 
good results in clusters but fails for absorption spec- 
tra of solids. Instead, more sophisticated approaches 
derived from Many-Body Perturbation Theory (MBPT) 
@j S B B have been able to reproduce, ab initio, 
the effect of the electron-hole interaction in extended sys- 
temSjiiot least thanks to an explicit long-range contribu- 
12l | . The latter strongly influences spectra like 



11 



tion 

optical absorption or energy loss, especially for relatively 
small momentum transfer. 

Here we show that this kernel is even able to reproduce 
the hydrogen-like excitonic scries in the photoemission 
gap of a rare gas solid. However the kernel has a strong 
spatial and frequency dependence, and its evaluation re- 
quires a significant amount of computer time. We there- 
fore tackle the question of a parameter-free, but quick 
TDDFT calculation of excitonic effects in solids, which 
has been so far an unsolved problem, and show that a 
much more efficient formulation can indeed be achieved. 
In particular we demonstrate how it is possible for a wide 
range of materials to obtain good absorption spectra in- 
cluding excitonic effects with a static kernel leading in 
principle to a Random Phase Approximation (RPA)-like 
scaling of the calculation with the number of atoms of 
the system. 

Atomic units are used throughout the paper. The vec- 
torial character of the quantities r, k, G, q (where k and 
q are vectors in the Brillouin zone, and G is a reciprocal 



lattice vector) is implicit. Only transitions of positive 
frequency (i.e. resonant contributions), which dominate 
absorption spectra, are considered throughout. 

Let us first concentrate on the absorption spectrum of 
solid argon. The low band dispersion, together with the 
small polarizability of the solid, conjures a picture where 
the electron-hole interaction is very strong and gives rise 
to a whole series of bound excitons below the intcrband 
threshold. As in the optical spectra of other rare gas 
solids, the first exciton is strongly bound (in argon by 
~ 2 eV), falling in the class of localized Frenkel [l6| exci- 
tons. Closer to the continuum onset at 14.2 eV, one finds 



more weakly bound Mott-Wannier [l7l | type excitons in 
a hydrogen-like series. In the ab initio framework, such a 
complex spectrum is typically described by the solution 
of the four-point (electron-hole) Bethe-Salpeter equation 



(BSE) 0, 14, 1^. In FigUwc show the optical spectrum 
of solid argon calculated within the BSE approach, and 
within TDDFT both using TDLDA ^ and the MBPT- 
derived kernel Q ■ The agreement of the BSE curve with 
experiment (line-circles) |2l| (and with previous BSE cal- 
culations [23) is good, concerning both position and rel- 
ative intensity of the first two peaks. It should be noted 
that the experiment shows double peaks due to spin-orbit 
splitting, which is not taken into account in our calcula- 
tions. The latter yields the singlet excitons that should 
essentially relate to the hole with j = 1/2 and be com- 
pared with the n' peaks. Besides the spin-orbit splitting, 
the pseudopotential approximation as well as the con- 
struction of a static W from LDA ingredients contribute 
to the remaining discrepancy with experiment. In spite of 
these limitations, the n' = 3 peak can also be detected, 
although the 2048 k-points used to calculate the spec- 
trum are not sufficient to discuss it quantitatively, nor to 
describe the higher peaks. Instead, the first two peaks 
require less k-points and, as can be seen in the inset, are 
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Figure 1: Absorption spectrum of solid Ar. The BSE (dashed) 
and TDDFT using kernel of Ref . (solid) are compared (only 
the n' singlet exciton series) with experiment [2l|. TDLDA 
is given by the points. Main panel: calculation with 2048 
k-points. Inset: 256 k-points 



already well reproduced with 256 k-points. In the follow- 
ing we therefore concentrate on these two structures and 
perform all calculations with 256 k-points. 

The BSE impressively improves upon the TDLDA 
(dotted), which shows a structure-less broad curve, 
clearly missing the bound excitons. Instead, the kernel of 
Ref. [7| (full curve in Fig. [1]) leads to the same accuracy 
as the BSE, both for the Frenkel exciton and for the fol- 
lowing structures. This demonstrates the potential of the 
method and shows that the MBPT-derived kernel can be 
used to quantitatively predict the absorption spectra of a 
wide range of materials, including the insulating rare-gas 
solids. 

However, the method is still computationally relatively 
heavy. Indeed in the MBPT-derived TDDFT approach, 
right as for the BSE two-particle Hamiltonian, one has 
to evaluate the matrix elements F^^f^ of the statically 
screened electron-hole Coulomb interaction W, 



2Tra dridr2^Uri,r2)W{ri,r2)^t'{ri,r2) (1) 



pBSE 



where the product ^t{ri,r2) = <j>vk{ri)(l)l^.^q{r2) of two 
KS wavefunctions is a generalized non local transition 
term; here t is an index of transition with momentum 
transfer q, i.e. t = {vckq\, from valence vk to conduction 
ck + q states, a = 2/{Nki^o) with Nk number of k- 
points and Oo volume of the unit cell. The calculation 



of F^^f^ scales with the number of atoms Nat as NrNf ~ 
N^t, where is the number of points in real space, and 
Nt is the total number of transitions. Following Ref. 
0, one then constructs an approximate kernel /^f'"^ = 
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-AEt) " {u; + ir]-AEt') 
where (ri ) = $4 (ri , ri ) and AEt are differences be- 
tween quasi-particle (QP) eigenvalues, since is an 
approximation to the "many-body" kernel that has 
to be used in conjunction with an independent particle 
response function xo built with QP energies instead of 
Kohn-Sham (KS) ones as in pure TDDFT. This kernel 
simulates hence to good approximation the electron-hole 
interaction that is described by the BSE Q. 

Even though this construction can be optimized Q 
the method is at least an order of magnitude slower than 
an RPA calculation. In the following we show how this 
problem can be overcome. 

We concentrate on the irreducible polarizability P that 
yields via the bare Coulomb interaction v the reducible 
polarizability x from the matrix equation x = P + P^x, 
and the inverse dielectric matrix e^^ = 1-1- vx- AH 
quantities are functions of q and frequency cj, and ma- 
trices in G,G' . Absorption spectra are then obtained 
from Abs(a;) = liin Im |l/eQQ^((7,a;)}. The polarizabil- 
ity P is determined from the screening equation P = 
Xo + XgI^qP- this equation we can now insert to the 
left and right of the identity 1 = XX'^ = X'^X, 
providing that X is a non-singular function. This yields 



P 



Xo + XoX-'TX-'P 



(2) 



where T = Xf"^^X. We choose a matrix of 
the form X = a (7t([jj)<I>t(r)<I>( (r'), where gt{(^) 
is an arbitrary function. The term T contains 
an explicit sum over matrix elements ^tddft _ 
Aira J dridr2^l{ri)f^^{ri,r2,LLi)^t{r2) in a basis of 
transitions $t, namely 

r(r, /, ^) = « ^ gt{u;)Mr)F^t'''''^gt' {u;)^ (r'). (3) 



The exact F™^^'^ is of course not known. However, in 
the spirit of Refs. 0, 0] we now replace the unknown 
matrix elements F™^^'^ with the BSE ones, given by 
Eq.([T]). With this mapping, T is approximated as 
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Mir')gt,{ij) = X^W^'X 



(4) 
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Figure 2: Main: illustration of the different choices for Eq. ((5| . 
In the inset: same choices for the absorption spectrum of Si. 



where we have defined a three-point right and left X op- 
erator as: X'^(ri]r2r2']uj) = a X)* 5t(^)*t(^i)^t (?'2, ^2') 
and ^X{rirv]r2]uj) = a X)* 5't('^)^t' ?'i')^t' (^2)- 
Here it is important to underline that i^'TDDFT ^-^^ 
structed as matrix elements of the local $t(r), whereas 
F^^^ are matrix element of the non-local ^t{r,r'). In 
fact the mapping ^ is not an exact operation, because 
Ff^f^ cannot be expressed as a matrix element (between 
$t and ) of a single f^^^ for all i , t' 0, 0, [H] . Therefore 
flf = x-^T^^X-^ can be different from f^^^, and the 
quality of the resulting spectra will depend on the choice 
of X. 

If a certain freedom in the choice of gt can be exploited, 
one may find approaches that boost the computational 
efficiency with respect to f^^'"^, i-e. the expression of 

In the following we will first illustrate, with the exam- 
ple of bulk Silicon and solid Argon, how different choices 
for gt{^) can lead to very similar spectra. 

We label with calligraphic letters the different choices 
A, B, C, V that stand for: 



A) gt{oj) = l/{uj~AEt+iv), (i.e X = xo) 

B) gt{uj) = lm{l/{uj~/^Et+ii^)} 

C) gt{oj)^l/AEt ; V) gt{oj) 



(5) 
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The first choice {A) defines nothing but the case X = 
XO; as proposed in Ref. Q and leading to /^c"^ above; 
in the second case (B) only the imaginary part is taken 
from the denominator of the independent particle polar- 
izability (very localized function in frequency); the cases 
(C) and (2?) describe simple static choices for gt{(^)- 

The inset of Figl2] shows the optical absorption of bulk 
Sihcon calculated with the BSE and within TDDFT, us- 
ing these mapping kernels {A,B,C,T>); the TDLDA re- 
sult is also shown in order to emphasize the little differen- 
cies among the mapping kernels, compared to the huge 
improvements of {A, BjC,!)) with respect to TDLDA. 
The description of the optical absorption of Argon is a 



much more stringent test. In Fig [2] we see that all the 
different kernels {A being slightly better than the oth- 
ers) arc able to well reproduce the excitonic series and to 
strongly improve upon the TDLDA result (dotted curve) . 
This is especially surprising for choices (C) and (V): 
bound excitons have up to now only been obtained using 
either the full, strongly frequency dependent kernel [9| 
or a frequency dependent long-range model {a + Pu'^) / q'^ 
24| , whereas a static scalar model can at the best yield 
one single bound exciton, with largely overestimated in- 
tensity, by tuning appropriately two model parameters 
[2^ . Our excellent results of FiglT] show, for the first 
time to the best of our knowledge, that even a static 
parameter-free two-point kernel is able to reproduce a 
scries of strongly bound excitons. 

It is now crucial to understand and hence predict 
the performance of the various choices, and to elucidate 
whether one can choose any possible X . To this aim we 
start from the four-point Bethe Salpeter equation *P 
Xo+ \oW *P, and contract the left and right indices. We 
obtain hence P ^ Xo + Xo^ "^-^i where we have defined a 
three-point right xo as Xo = S:o{ri,ri;r2,r2';Lu) and a 
three point left polarizability = ^P(ri, ri*; r2, r2; w). 
Now, inserting the identity 1 — XoXo ^- we obtain: 



P 



Xo 



XoXo ^xlW 



'P. 



(6) 



On the other hand using the mapping Q in eq.Q, we 
obtain the approximate polarizability 



P' 



.cff 



1 noff 



XQ + XoX~'X'W'XX~'P 



(7) 



If TDDFT is to reproduce the BSE results, P"^ result- 
ing from ^ must be equal to P of (l6|). 

It should be noted that in principle the matrix X can 
be chosen differently for the left and for the right side of 
W in eq. ([7]) and eq.®. Concentrating first on the left 
side, the choice X = xo, i-C 5t(w) = l/(^ + irj ~ AEt) 
recovers exactly the left side of W in eq. ([6]). 
The right side is still to be optimised. The comparison 
between dH) and ([7]) suggests to choose X — P. Of course 
this is not the solution of the problem, since i) P is the 
quantity we are looking for and ii) P cannot be expressed 
as a sum over KS transitions respecting the ansatz for X. 
Hence, one can only try to find a good guess. Again, A, 
with X = Xo seems a good choice. In fact, in a solid 
the joint density of states calculated in GW is very close 
to the density of transition energies calculated from the 
BSE; i.e. xo from GW and P from the BSE have a very 
similar distribution of poles (2^. Concerning the other 
choices, it is useful to note that P^^^ and P have the same 
poles; the same statement holds for X*-^-* and X. If, in 
Eq.© the poles of X^'^^ cancelled with the zeroes of X~^, 
and no new poles were introduced, one would just find the 
poles of P°^ in the right side of W in Eq. ([7]) , right as for P 
in ([6]). However, X^^ has poles that lie between the poles 
of X . These new poles are not problematic for energies in 
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the continuum, but they can lead to spurious structures 
when they appear isolated, i.e. in the bandgap. It turns 
out that this effect is particularly strong when the poles 
of X are in the vicinity of the bound excitons. This is 
for example the case when one chooses X = Xo^ 
gtioj) = l/{uj + i?7 - AE'f •^), AiJf being the difference 
between KS eigenvalues): indeed, the pink circles in Fig. 
[2] show the bad performance of that choice.. 

We have, in fact, verified that the spectra are gener- 
ally very stable as long as we choose an X that (i) either 
does not have any poles (static choices); or (ii) has poles 
in the continuum (like xo); or (iii) has poles at very low 
energies, much lower than all poles of P®*^. This con- 
firms that a wide range of choices for X is indeed possi- 



ble. Moreover, this observation is valid for a wide range 
of materials: we have performed the same test calcula- 
tions for the prototype materials diamond and SiC, with 
similar conclusions. 

As pointed out above, the aim is to avoid the unfavor- 
able scaling of the calculations, determined essentially by 
the evaluation of F^^^^ via Eq. H]) . Choice (2?) is of course 
particularly simple and promising. In fact even when it is 
used as it is in (jH) , the static choice V leads to a speedup 
with respect to choice A 0] ■ More importantly, it allows 
one to recombine the sums and integrals in Eq.Q in a 
more convenient way. The latter equation, once {V) is 
chosen, can in fact be written as 



r°*^(<7,G,G') = -47Ta^Y.^o,o'{q)fdrdr'dfdr'e-''' '^e'''' '^'Ak{r, r)Bk-q{r', r)Ak+q{f, r')Bk+q 

kq GG ^ 



(8) 



where Ak{r,r') = Y.vKk{T)uvk{r'), Bk{r,r') = 
Sc'"cfc(''')'"cfe(^')5 with the u's representing the periodic 
part of the KS wavefunction (l)vk{r) = e~'^^'^Uyk{i')- 
Wg 5'(g) is the reciprocal space Fourier transform of the 
statically screened Coulomb interaction, with q the dif- 
ference between two k-points in the Brillouin zone. For 
q we have the special case of vanishing momentum 
transfer q (e.g. for optical absorption); ([8]) is the gen- 
eral formula valid for any q in order to treat also, e.g., 
electron energy loss or inelastic X-ray scattering. 

The scaling of Eq.® is in principle iV^^, but with a 
clearly dominant contribution given by the spatial inte- 
grals, which scales as N'^\n{Nr) < N^t I!!- Note that 
iV^t is the scaling of the construction of xo itself [23|. 
In other words, this formulation offers the possibility to 
determine absorption spectra including excitonic effects 
with a workload comparable to the RPA. 

In conclusion, we have calculated the absorption spec- 
tra of solid argon, both by solving the Bethe-Salpeter 
equation and by time-dependent density functional the- 
ory using a MBPT-derived mapping kernel. Both meth- 
ods yield results in good agreement with experiment and 
reproduce well positions and relative intensities of the 
peaks, with a drastic improvement over TDLDA results. 
We have then introduced a method that allows one to de- 
rive a variety of approximations for the TDDFT kernel; 
these can be used to tune computational efficiency while 
maintaining most of the precision of the original formula- 
tion. The method has been tested for solid argon, silicon, 
diamond and silicon carbide. The good results, in turn, 
have allowed us to propose a reformulation of the ker- 
nel (Eq.®) that leads to a TDDFT calculation with the 
same quality of the BSE, but with an RPA-likc scaling 



< N'^^, rather than N^^. This, we believe, can consti- 
tute a real breakthrough for practical applications where 
a low computational effort - that characterizes TDDFT 
- and a precise description of many-body effects - like in 
the BSE - are required. 
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